/**
 * 
 */
package ikelib;

/**
 * 加速度の時系列データから、速度、変位に変換する。
 * 
 * 基線の補正、最大加速度の補正を行う。
 * 
 * 
 * @author Toshiharu Ikeda <ikeda.104@gmail.com>
 *
 */
public class IntegrateAccelaration extends Integration {
	/**
	 * 時刻のステップ　[sec]
	 */
	private double[] Acs = null;	// 補正された加速度
	private double[] Vel = null;	// 補正された速度
	private double[] Dsp = null;	// 補正された変位
	private double T = 0.0; 	// 積分区間
	private double a0 = 0.0; // 基線（base line）補正係数 a0
	private double a1 = 0.0; // 基線（base line）補正係数 a1
	
	
	
	public double[] getAcs() {
		return Acs;
	}


	public double[] getVel() {
		return Vel;
	}

	public double[] getDsp() {
		return Dsp;
	}

	public static void main(String[] arg){
		double[] test = new double[100];
		for(int i=0;i<100;i++){
			//test[i] = Math.sin(2*Math.PI * i/100);
			if( i>20 && i<40) test[i] = 1;
			else if(i>50 && i<60) test[i] = -2;
			else test[i] = 0;
		}
		IntegrateAccelaration intd=new IntegrateAccelaration(0.1);
		intd.Calc(test);
		System.out.println("end");
	}
	
	/**
	 * コンストラクタ
	 * @param timestep データの時刻幅 [sec]
	 */
	public IntegrateAccelaration(double timestep){
		super(timestep);
	}
	
	/**
	 * 加速度データを与えて計算
	 */
	@Override
	public boolean Calc(double[] acs){
		return Calc(acs,0.0);
	}
	/**
	 * 加速度データ,オフセットを与えて計算
	 * @param acs
	 * @return
	 */
	@Override
	public boolean Calc(double[] acs,double ofs){
		// ＤＣオフセットをキャンセルする
		Offset(acs, ofs);
		Statistics acsSta = CalcStatistics(acs);
		// 最大値、最小値、などを保存
		// Ｔを計算
		T = TimeStep*acsSta.Length; 
		// 速度と変位を数値積分
		//   a１の積分も同時に計算する。
		int len = acsSta.Length;
		double[] itg = new double[len];
		Vel = new double[len];
		Dsp = new double[len];
		Acs = new double[len];
		for(int i=0;i<len;i++){
			if(i==0){
				Vel[i] = 0.0;
				Dsp[i] = 0.0;
				itg[i] = 0.0;
			}
			else{
				// 速度を計算
				Vel[i] = Vel[i-1] + (acs[i-1] + acs[i])*TimeStep/2.0;
				//　変位を計算
				Dsp[i] = Dsp[i-1] + Vel[i-1]*TimeStep + (acs[i-1]/3 + acs[i]/6)*TimeStep*TimeStep;
				// integpartを求める
				double t = TimeStep*i;
				itg[i] = Dsp[i]*(3.0*T*t*t- 2.0*t*t*t);
			}
			//System.out.print(String.format("%f\t%f\t%f\t\n",acs[i],Vel[i],Dsp[i]));
		}
		//System.out.println("Hosei");
		// a1を求めるための台形則の積分を実行する。
		double a1integpart = TrapezoidalIntegration(itg,TimeStep);
		
		// 基線補正用係数
		a1 = 28.0/13.0/(T*T)*(2.0*Vel[len-1] - 15/(T*T*T*T*T)*a1integpart);
		a0 = (Vel[len-1] - a1*T*T/2.0)/T;
		// 基線補正
		for(int i=0;i<len;i++){
			double t = TimeStep*i;
			Acs[i] = acs[i] - (a0 + a1*t);
			Vel[i] = Vel[i] - (a0*t + a1*t*t/2.0);
			Dsp[i] = Dsp[i] - (a0*t*t/2.0 + a1*t*t*t/6);
			//System.out.print(String.format("%f\t%f\t%f\t\n",Acs[i],Vel[i],Dsp[i]));
		}
		
		// 最大振幅が一致するように補正する　→最大値が一致するように補正する。
//		double amp = acsSta.Max-acsSta.Min;
//		Statistics acsSta2 = CalcStatistics(Acs);
//		double amp2 = acsSta2.Max - acsSta2.Min;
//		Acs = Amp(Acs,amp/amp2);
		double amp;
		if(Math.abs(acsSta.Max) > Math.abs(acsSta.Min)){
			amp = Math.abs(acsSta.Max);
		}
		else{
			amp = Math.abs(acsSta.Min);
		}
		Statistics acsSta2 = CalcStatistics(Acs);
		double amp2;
		if(Math.abs(acsSta2.Max) > Math.abs(acsSta2.Min)){
			amp2 = Math.abs(acsSta2.Max);
		}
		else{
			amp2 = Math.abs(acsSta2.Min);
		}
		Acs = Amp(Acs,amp/amp2);
		
		return true;
	}
	
	
	
}
